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Abstract. In this paper we develop a finite-difference scheme to approximate 
radially symmetric solutions of the initial-value problem with smooth initial 
conditions 

— - v 2 «, - 0— (V 2 ™) + 7— + m 2 w + G (w) = 

subject to : f w (x,0) = <p(x), iEl) 

dw 

dt 



^ ' ^(2,0) = tp(x), xeD 



in an open sphere D around the origin, where the internal and external damp- 
ing coefficients — and 7, respectively — are constant, and the nonlinear term 
has the form G'(w) = w p , with p > 1 an odd number. The functions <p and tp 
are radially symmetric in D, and <f>, ip, rep and rip are assumed to be small at 
infinity. We prove that our scheme is consistent order 0(At 2 ) + 0(Ar 2 ) for G 1 
identically equal to zero, and provide a necessary condition for it to be stable 
order n. Part of our study will be devoted to compare the physical effects of 
and 7. 



1. Introduction 

Klcin-Gordon-like equations appear in several branches of modern physics. A 
modified sine-Gordon equation appears for instance in the study of long Josephson 
junctions between superconductors when dissipative effects are taken into account 
[I] . A similar partial differential equation with different nonlinear term appears in 
the study of fluxons in Josephson tramsmission lines [5] . A modified Klein-Gordon 
equation appears in the statistical mechanics of nonlinear coherent structures such 
as solitary waves in the form of a Langevin equation (see [3] pp. 298-309); here 
no internal damping coefficient appears, though. Finally, our differential equation 
describes the motion of a damped string in a non-Hookean medium. 

The classical (I + I)-dimensional linear Klein-Gordon equation has an exact 
soliton-like solution in the form of a traveling wave 0]. Some results concerning 
the analytic behavior of solutions of nonlinear Klein-Gordon equations have been 
established [SJ 111 [7]; however, no exact method of solution is known for arbitrary 
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initial- value problems involving this equation. From that point of view it is im- 
portant to investigate numerical techniques to describe the evolution of radially 
symmetric solutions of ([1]). 

It is worth mentioning that some numerical research has been done in this direc- 
tion. Strauss and Vazquez [5] developed a finite-difference scheme to approximate 
radially symmetric solutions of the nonlinear Klein-Gordon equation for the same 
nonlinear term we study in this paper; one of the most important features of their 
numerical method was that the discrete energy associated with the differential equa- 
tion is conserved. The numerical study of the sine-Gordon model that describes the 
Josephson tunnel junctions has been undertaken by Lomdahl et al. [5] . Numerical 
simulations have also been performed to solve the (1 + l)-dimensional Langevin 
equation [9]. 

In this paper we extend Strauss and Vazquez's technique to include the effects 
of both internal and external damping, and validate our results against those in [8]. 
Section [2] is devoted to setting up the finite-difference scheme; the energy analysis 
of our problem is also carried out. Numerical results are presented in Section [31 
followed by a brief discussion. 



2. Analysis 

Analytical results. The following is the major theoretic result we will use in our 
investigation. Here M(t) represents the amplitude of a solution of (p} at time t, 
that is 

M(t) = max \ w(x, t)\. 

X 

Theorem 1. Let j3 and 7 be both equal to zero, and let G'(w) = |u>| p ~ 1 u'. Suppose 
that <f> and ip are smooth and small at infinity. Then 

(1) If p < 5, a unique smooth solution of |T]) exists with amplitude bounded at 
all time [B] . 

(2) If p > 5, a weak solution exists for all time [10) . 

(3) For p > 8/3 and for solutions of bounded amplitude, there is a scattering 
theory; in particular, they decay uniformly as fast as M(t) < c(l + |i|)~ 3 ^ 2 

M- □ 

Finite-difference scheme. Throughout this section we will assume that the func- 
tions <f>{x) and ip(x) are smooth, of compact support, radially symmetric in the open 
sphere D with center in the origin and radius L, and that <f>, ip, r<f> and rip are small 
at infinity in D. Moreover, we will suppose that w(x,t) is a radially symmetric 
solution of (JTJ) . 

Let r = ||x|| be the Euclidean norm of x and let G'(w) — w p , for p > 1 an odd 
number. Setting v(r, t) = rw(r, t) for every < r < L and t G R, it is evident that 
v must satisfy the mixed- value problem 

d 2 v d 2 v dv n d 3 v 2 , . 

w - + ^t - ^ + m v + rG {v/r) = 



(2) subject to 



v(r, 0) = r<j)(r), < r < L 

— {r,0)=rijj(r), < r < L 
at 

v(0,t) = 0, t>0. 
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Proceeding now to discretize our problem, let a < L be a positive number with 
the property that (j) an d ip vanish outside of the sphere with center in the origin 
and radius a — e, for some e > 0. Let = r < r\ < ■ ■ ■ < Tm — a an d = 
to < t\ < ■ ■ ■ < tjv = T be partitions of [0,a] and [0,T], respectively, into M and 
N subintervals of lengths Ar = a/M and At = T/N, respectively. Denote the 
approximate value of v(rj,t n ) by vj. The finite-difference scheme associated with 
© is 

v n + l _ 2v n + p n-l ^ _ 2y n + ^ n + l _ g-1 

(At) 2 (Ar) 2 7 2At 

n+l q n+l , n+l\ / n— 1 q n— 1 j_ 



(3) /9-^± 2 M±i * + 

2At(Ar) 2 

where G(¥) = 1). 

Computationally, our method requires an application of Newton's method for 
systems of nonlinear equations along with Crout's reduction technique for tridiog- 
onal linear systems. 

An interesting property of this finite-difference scheme is that, for initial ap- 
proximations {wj} and {wj} with zero centered-difference first spatial derivatives 
at the origin, the successive approximations provided by the method have likewise 
centered-difference first spatial derivative equal to zero at the origin. This claim 
follows by induction using the facts that 

* y (Vl +1 + A*) - (vt 1 + A 1 ) 

that Wi = w n L 1 iff + v™i — 0, and the substitution = jw™ Ar. The induction 
hypothesis implies that v" + v™ 1 = for every n < k, whence the claim follows. 

Our last statement implies that for a smooth initial profile at the origin, the 
subsequent approximations yielded by our method will be likewise smooth. As a 
test case, it is worthwhile mentioning that we have successfully obtained numerical 
results to verify this claim using a Gaussian initial profile centered at the origin. 

Stability analysis. It is clear that © is consistent order 0(At 2 )+0(Ar 2 ) with © 
whenever G' is identically equal to zero. Moreover, in order for the finite-difference 
scheme to be stable order n it is necessary that 

At\ 2 At „ At 2 (At) 2 

A r ) 4 (Ar) 2 4 

To verify this claim, notice first that ([3]) can be rewritten as 

v] +1 - 2vj + v™- 1 s 2 vj v r ; +1 - v r ;- Y 



(At) 2 (Ar) 2 + 7 " 2At 

5^-5^ , m 2 L , n+1 , 



K +1 



2At(Ar) 2 2 1 J 

Define R = At/Ar. Let Vfi +1 = and V£ +1 = vj for each j = 0, 1, . . . , M 
and n = 0, 1, . . . , N - 1. For every j = 0, 1, . . . , M and n = 1, 2, . . . , N let VJ 1 be 
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the column vector whose components are V{j and VJj • Our problem can be written 
then in matricial form as 



where 



k \^ n+1 / 2 + RHl -h . 

I) 1 / I 1 I j ' 



At /3At<« 2 (At) 2 

1 + 7—-- T+^-^TT^ and 

2 2 (Ar) 2 2 



* = 1-7^ + ^W^. 

2 2 (Ar) 2 2 

Denoting the Fourier transform of each VP by VP, we obtain that 



(l - 2R 2 sin 2 | 



VP+ 1 = ( fe(i) V A "' 2 7 HQ I 



where 



At BAt . 2 £ 2 (^) 2 

h2 7 sm — h m 

2 (Ar) 2 2 2 



fc(f) = 1 + 7— + 2^ — 2 sin 2 | + m 2 -^^- and 



At 8 At 2 f , (At)" 
^ = i-T^T- 2 — — o-sin 2 J +m 2 ^-^- 
2 (Ar) 2 2 2 



The matrix A(£) multiplying VJ 1 in the last vector equation is the amplification 
matrix of the problem, which has eigenvalues given by 



1 - 2i? 2 sin 2 | ± x j ( 1 - 2i? 2 sin 2 | ) - h(£)k(Q 
A± = 



In particular, for £ = it the eigenvalues of A are 



1 - 2i? 2 ± J(l - 2i? 2 ) 2 - /i(7r)jfe(7r) 

A± = * . 

&(tt) 

Suppose for a moment that 1 — 2i? 2 < —k(Tr). If the radical in the expression for 
the eigenvalues of A(tt) is a pure real number then |A_| > 1. So for every n € N, 
1 1 -A™ 1 1 > I A- 1" grows faster than + nK2 for any constants i^i and K?,- A similar 
situation happens when the radical is a pure imaginary number, except that in this 
case | • | represents the usual Euclidean norm in the field of complex numbers. 

Summarizing, if 1 — 2i? 2 < — k(ir) then scheme ([3]) is unstable. Therefore in order 
for our numeric method to be stable order n it is necessary that 1 — 2i? 2 > —k(ir), 
which is what we needed to establish. 

Energy analysis. Assume that G: R — > E is continuously differentiable, and that 
w(x,t) is a solution of ([TJ in a domain D of R 3 . Moreover, we assume that Vw • n 
is zero near the boundary of D at all time, where n denotes the unit vector normal 
to the boundary of D. The Lagrangian associated with our nonlinear modified 
Klein- Gordon equation is given by 

£ = -{[ -f ) -\Vw\ 2 -m 2 w 2 \-G(w). 
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0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



Figure 1. Approximate radial solutions of JT} at successive times for (5 = 
and values of 7 = (solid), 7 = 5 (dashed) and 7 = 10 (dotted), G'(u) = u 7 , 
and initial data (j>(r) = h(r), ^(r) = h'(r) + h(r)/r. 



It is easy to derive the following expression for the total energy associated with our 
nonlinear dissipative Klcin-Gordon-like equation: 

(4) E(t)= [[[ (if^V + W + ^+GHL,. 



D 



2\dt 2 1 ' 2 
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Figure 2. Approximate radial solutions of JT} with G'(u) at t = 0.2, for initial 
data tj>(r) = and ip(r) = W0h(r), /3 = and 7 = (solid), 7 = 5 (dashed) 
and 7 = 10 (dotted). 



Proposition 2. The instantaneous rate of change with respect to time of the total 
energy associated with the PDE in ([!} is given by 
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Figure 3. Approximate radial solutions at successive times of the undamped 
(solid) and the damped nonlinear Klein-Gordon equation with 7 = and 
P = (dashed), p = 0.0001 (dash-dotted) and P = 0.0002, nonlinear term 
G 1 (u) = u 7 , and initial data <j>(r) = h(r) and ip{r) = h'(r) + h(r)/r. 



Proof. Taking derivative on both sides of Equation (QJ , we obtain that 
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G'(u) = u 3 G'{u) = u 3 




-7 5 

1 1 1 i 1 1 1 1 1 1 _gl 1 1 1 1 1 1 1 

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



Figure 4. Approximate radial solutions of JT} with G'(u) at t = 0.2, for initial 
data </>(r) = 0, ip(r) = 100/l(r) (left column) and (p(r) = h(r), tp(r) = (right 
column), 7 = and f3 = (solid), = 0.0001 (dashed) and /3 = 0.0002 
(dotted). 



On the other hand, from Green's first identity we see that 



dw. 



D 



at 



dw 



Of 



dw d 



— dx= __( Vw .n) da- 



<>D 



dt dt 



D 



VI — 

dt 



dx. 



The surface integrals in these last two equations are equal to zero, whence the result 
follows. □ 
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0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



G' (u) — u 9 G' (u) — sinh(5tt) — bu 




Figure 5. Approximate radial solutions of JTJ at t = 0.2, for 7 = 5, initial 
data <t>(r) = and tp(r) = 100fe(r), and f3 = (solid), = 0.0005 (dashed) 
and p = 0.005 (dotted). 

It is worthwhile noticing that if /3 and 7 are positive then the total energy is 
decreasing in time. Also, if /3 and 7 are both equal to zero then the energy is 
conserved. Finally, if j3 is zero then the expression of E'(t) coincides with the one 
derived in [12] . 
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Figure 6. Approximate value of solutions to JTJ near the origin vs. time for 
different nonlinear terms, and initial conditions <p(r) = and ip(r) = 100/i(r). 
Left column: f) = and 7 = (solid), 7 = 10 (dashed) and 7 = 20 (dotted); 
right column: 7 = and = (solid), = 0.001 (dashed), = 0.0025 
(dashed-dotted) and = 0.005 (dotted). 



Let us assume now that G is nonnegative. The total energy in this case is likewise 
nonnegative and the integral of every term in (QJ is bounded by y / 2E(t)/m. In 
particular, this last statement implies that the integral of w 2 at time t is bounded 
by E(t). For those times t for which E(t) is finite (and particularly for the case 
when /? and 7 are both equal to zero), this means that w is a square- integrable 
function in the first variable at time t. 
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G'(u) = u 3 G'(u) = u 3 




Figure 7. Total energy vs. time for initial data <f>(r) = and ip{r) = 100fo(r). 
Left column: = and 7=1 (solid), 5 (dashed), 10 (dotted). Right column: 
7 = and /3 = 0.0005 (solid), 0.001 (dashed), 0.005 (dotted). 

Let G'(w) — w p with p > 1. Assuming that w is a radially symmetric solution 
of the damped nonlinear Klein-Gordon equation in a sphere D with center in the 
origin and radius L, and using the transformation v(r,t) = rw(r,t) the energy 
expression adopts the form E(t) — 47r_E (i), with 

< 5 » E ^[[M + \^ + ^ 2+r '-' a ^ 
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Table 1. Relative differences of externally damped solutions to ((TJ with re- 
spect to the corresponding undamped solution at different time steps. 



Time step 




Relative differences 




n 


7 = 0.1 


7 = 0.5 


7=1 


7 = 5 


7 = 10 





0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


20 


0.0028 


0.0142 


0.0283 


0.1395 


0.2693 


40 


0.0103 


0.0509 


0.1006 


0.4491 


0.7706 


GO 


0.0167 


0.0821 


0.1611 


0.6579 


0.9573 


80 


0.0192 


0.0942 


0.1836 


0.6954 


0.9387 


100 


0.0200 


0.0977 


0.1896 


0.6994 


0.9308 



The instantaneous rate of change of energy is given by E'(t) — AirE' (t), where 
f L \ a fd 2 v ldv\ 2 



\ drdt r dt J \dt 




It is possible to reproduce now the argument in [8] to show that for every t and 
nonzero r, \w(r,t)\ < v / 2E (t)/r. This means in particular that if a solution were 
unbounded, it would have to be unbounded at the origin. 

The discrete energy is given by 
Eq 

Ar 2 ^ \ Al j 2 ^ \ Ar J \ Ar 




E 



t 2(jAr)p- 

This expression is obviously consistent with ([5|). Moreover, taking the difference 
between Eq / Ar and E^ 1 / Ar and simplifying after using ©, it can be shown that 

EE K- 1 _ B V ( V J__Z^1 \ A< +1 -^n-(^ + i 

At f^i \ 2At ) \ At(Ar) 2 





2At 



Ar. 



For f3 = this expression provides us with a consistent approximation to the 
instantaneous rate of change of energy. Numerical results demonstrate that en- 
ergy decreases as a function of time for /? > 0, which is in agreement with the 
corresponding instantaneous change of energy as a function of time. 

3. Numerical results 

The numerical results presented in this section correspond to approximate solu- 
tions of the dissipative, nonlinear, modified Klein-Gordon equation 

gH_V^-^(V^ )+7 g + „ + G »=0, 
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Table 2. Table of relative differences of externally damped solutions of JTJ 
with respect to the corresponding undamped solution at t = 0.2. 



Nonlinear Term 




Relative differences 




G'(u) 


7 = 0.1 


7 = 0.5 


7 = 1 


7 = 5 


7 = 10 





0.0098 


0.0478 


0.0923 


0.3642 


0.5631 


u 3 


0.0097 


0.0477 


0.0929 


0.3528 


0.5554 


u 5 


0.0137 


0.0665 


0.1287 


0.4024 


0.6418 


u 7 


0.0171 


0.0833 


0.1618 


0.5068 


0.7819 


q 

vr 


0.0204 


0.0999 


0.1728 


0.5736 


0.8488 


sinh(5«) — 5u 


0.0263 


0.1377 


0.2518 


0.6284 


0.8813 



Table 3. Relative differences of internally damped solutions to |(TJ with re- 
spect to the corresponding undamped solution at different time steps. 



Time step Relative differences 



n 


p = io- b 


p = io-'° 


p = 10~ 4 


P = 0.0005 


p = 0.001 





0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


20 


0.0005 


0.0054 


0.0517 


0.2188 


0.3682 


40 


0.0030 


0.0300 


0.2640 


0.8281 


1.1457 


60 


0.0156 


0.0996 


0.1493 


1.6536 


1.1460 


80 


0.0102 


0.0970 


0.7242 


1.1530 


1.2138 


100 


0.0080 


0.0772 


0.5751 


1.0406 


1.1435 


Table 


4. Table of relative differences 


of internally damped solutions of JTJ 


with respect to the corresponding undamped solution at t = 0.2. 




Nonlinear Term 


Relative differences 




G'(u) 


P = 10- 


b p = io- 


p = 10~ 4 


P = 0.0005 


p = 0.001 





0.0003 


0.0027 


0.0242 


0.0859 


0.1326 


u 3 


0.0003 


0.0032 


0.0289 


0.1040 


0.1621 


u 5 


0.0011 


0.0105 


0.0948 


0.3374 


0.5043 


u 7 


0.0023 


0.0224 


0.1825 


0.5663 


0.7327 


u 9 


0.0041 


0.0397 


0.3133 


0.7318 


0.9256 


sinh(5w) — 


5u 0.0063 


0.0577 


0.4717 


0.9403 


1.1007 



obtained using a tolerance of 10 5 and a maximum number of 20 iterations on 
every application of Newton's method. The space and time steps are always fixed 
as Ar = At = 0.002. 

External damping. Throughout this section we fix P — 0. 

Let us start considering the problem of approximating radially symmetric solu- 
tions of ([T]) with G'(u) — u 7 , and initial data <p(r) — h(r) and ip(r) = h'(r) + h(r)/r, 
where 

1 



h(r) 



5 exp <j 100 
0, 



1 - 



1 - (lOr- l) 2 



if < r < 0.2, 
if 0.2 < r < 0.4. 



We have plotted numerical solutions of this problem for several values of 7. The 
graphical results are presented in Figure Q] for 7 = 0,5,10. We observe first of all 
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that the solutions of the damped nonlinear Klein-Gordon-like equation correspond- 
ing to small values of 7 are consistently similar to those of the undamped case. To 
verify this claim quantitatively, we consider the approximations Vq and C™ to the 
undamped and damped cases, respectively, and compute the relative difference in 
the ^.Aariiorm via 

1 1 ^-.n -n I 
\\V \ \2,Ax 

(here we follow [T3]). The relative differences for several values of 7 at consecutive 
time steps are shown in Table [TJ We observe that the difference between the 
solutions of the nonlinear Klein-Gordon-like equation with damping coefficient 7 
and the corresponding undamped equation can be made arbitrarily small by taking 
7 sufficiently close to 0. 

We wish to corroborate this pattern for different nonlinear terms and a different 
set of initial conditions. With this objective in mind, Figure [5] depicts numerical 
solutions of ([T]) with 7 = 0,5,10, nonlinear terms G'(u) = 0, u 3 , u 5 , u 7 , u 9 , and 
sinh(5w) — 5w, initial conditions <j){r) = and ip( r ) — 100/i(r), and values of 7 = 
0, 5, 10. More accurately, Table [2] provides relative differences of these solutions 
at t = 0.2 (for the nonlinear functions listed above and varying values of 7) with 
respect to the corresponding undamped solution, for a wider selection of values of 
the parameter 7. 

It must be mentioned that, as it was expected, the total energy was invariably 
decreasing for positive values of 7, and increasing for negative values. For the 
value 7 = 0, the rate of change of energy is equal to zero and our numeric results 
agree with [8] . Experimental results show that small values of 7 correspond with 
small values of the discrete rate of change of the energy. This last observation 
corroborates stability of our method. 

We also observe that the amplitude of solutions corresponding to positive val- 
ues of 7 tend to decrease as time or 7 increases. Figure [2] partially corroborates 
that behavior. We have computed solutions corresponding to negative values of 7 
(graphs not included in this paper) and have verified that the amplitude of solutions 
increases with time and with \<y\. 

Finally, we have obtained graphs of the energy Eq vs. time for G'{u) = u 3 , u 5 , u 7 , 
initial data 0(r) = and ^(r) = 100/i(r), and values of 7 = 1,5,10. The results 
(depicted in the left column of Figure[7]) show a loss in the total energy as a function 
of time. 

Internal damping. Consider first the case when 7 equals zero. Figure [3] shows nu- 
merical solutions of ([I} at consecutive times, for initial data (j)(r) — h(r) and ip(r) = 
h'(r) + h(r)/r, with nonlinear term G'(u) = u 7 , and values of /? = 0,0.001,0.003. 
We observe that small values of f3 produce results similar to those of the corre- 
sponding undamped case. To corroborate this claim, we appeal once more to the 
relative differences in the ^2,A:r" nornl of dissipative solutions with respect to the 
non-dissipative one. The results are shown in Table [3] The results evidence the 
continuity of solutions with respect to the parameter /3 for this particular choice 
of nonlinearity, providing thus numerical support in favor of the stability of our 
method. 

We want to establish now the continuity of our method for several nonlinear 
terms at a fixed large time. In order to do it, Figure |4] shows the numerical solutions 
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of ([I]) at time t = 0.2, for the nonlinear terms G'(u) — u 3 ,u 5 ,u 7 , for two different 
sets of initial conditions: <f>{r) = and tp( r ) — 100(h), and <f)(r) — h(r) and 
ip{r) = 0, and values of = 0,0.0001,0.0002. The graphs in this figure, together 
with the analysis of relative differences in the £2, Ax-norm supplied in Table [4] for 
the first set of initial conditions, evidence the continuity of the numerical solution 
given by our method with respect to the parameter /3 for different nonlinearities. 

We now consider the case when 7 is nonzero. We use G'(u) = 0, u 3 , u 5 , u 7 , u 9 , 
and sinh(5w) — 5u, initial data <f>(r) = and ip( r ) = 100/i(r), and time t = 0.2. 
Figure [5] shows numerical solutions of (p} for values of /3 = 0,0.0005,0.005. The 
solutions for smaller nonzero values of j3 are indeed closer to the corresponding 
internally undamped solution, while the larger values of (3 spread out the internally 
undamped solution at the same time that the maximum amplitude is decreased. 

In order to study the time behavior of the solutions near the origin we have 
included Figure [SJ using initial data <f>(r) = and ip(r) = 100/i(r), the nonlinear 
terms G'(u) — u 3 , u 5 and u 7 , different choices of values for /3 and 7, and Ar = At = 
0.002. The left column shows the time-dependence of solutions at the origin for 
j3 = and three positive values of 7, whereas the right column shows similar results 
for 7 = and three positive values of /3. We observe that the value of solutions at 
the origin for large times is always approximately equal to zero for /3 = 0, which is 
in agreement with our experience of the (1 + l)-dimensional case. 

Finally, Figure [7] shows the graphs of the energy Eq vs. time for G'(u) = 
u 3 ,u 5 ,u 7 , initial data <j)(r) =0andV(r) = lOOh(r), and values of /3 = 0.0005,0.001,0.005. 
The results (depicted in the right column) show a loss in the total energy as a func- 
tion of time. It is clear that the rate at which the total energy is lost due to internal 
damping is greater than the corresponding rate due to external damping. 

4. Discussion 

A numerical method has been designed to approximate radially symmetric so- 
lutions of some dissipative, nonlinear, modified Klein-Gordon equations with con- 
stant internal and external damping coefficients (3 and 7, respectively. Our finite- 
difference scheme is in general agreement with the non-dissipative results presented 
in [8]. The method is consistent 0(At 2 ) + 0(Ar 2 ), conditionally stable, and contin- 
uous with respect to the parameters /3 and 7; as expected, the total energy decays 
in time for positive choices of the parameters. The corresponding scheme to ap- 
proximate the total energy of the system is consistent and has the property that 
the discrete rate of change of the discrete energy with respect to time approximates 
the corresponding continuous rate of change for /3 = 0. 

Several conclusions can be drawn from our numerical computations. First of 
all, we have seen that both internal and external damping tend to decrease the 
magnitude of solutions, as it was expected. Our results clearly exhibit the dispersive 
effects of the parameter (3 and the dissipative effects of 7. Our energy computations 
evidence the fact that the rate at which the energy is dissipated by the internal 
damping is faster than the corresponding rate of external damping. Finally, we 
observe that the effect of the nonlinear term in the temporal behavior near the 
origin is to increase the number of oscillations as the degree of the nonlinearity 
is increased. Invariably, the solutions of the dissipative modified Klein-Gordon 
equation converge in time to the trivial solution. 
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